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We designed and implemented our own versatile simulation software in order to understand the 
velocity changes and track patterns observed in slow moving vortex lattices. The data was obtained 
from time series of STM images on NbSe2 in a magnetic field range of 250 — 750 mT. The main thrust 
is to explore possible driving mechanism and to test the effect of various defect configuration on 
vortex velocity and tracks. The simulation uses the full vortex-vortex and vortex-defect interaction. 
It allows for periodic boundary conditions as well as a repulsive sample edge. In early versions the 
time intervals for recalculating the force were estimated for individual vortices. This, however, led 
to a positional 'noise' of unacceptable temperature of the vortex lattice and was hence replaced by a 
global time step control. We were able to produce similar track and velocity patterns as well as local 
lattice distortions near point defects as observed by STM. A detailed analysis is, however, beyond 
the scope of this paper. 



INTRODUCTION 



II. THE SIMULATION 



The dynamic behavior of vortices in type-II- 
superconductors has been extensively studied both 
in experimen1pii2ii^i^i^i^ and simulation. The simula- 
tions usually fall in one of two categories: Ginzburg- 
Landau-TheoryZi^i^ii^iiiii^ii^ or a 2D-molecular dynam- 
ics modeli^ii^ii^iii. The former is used to determine de- 
tails of the vortex configuration, usually for a small area 
of the sample, from first principles. The second is used 
to study the behavior of a large number of vortices. The 
work presented in this paper falls within the second cat- 
egory. Among the observed effects which are sought to 
be explained are the peak effect of the critical current 
in current driven vortex lattices as well as the dynamic 
behavior of the vortex lattice as imaged by several tech- 
niques such as magnetic decoration techniques'^, Lorentz 
microscopy^^ or scanning tunneling microscopy-'-^. 



The main purpose of the simulations described in this 
paper is to compare the results to the behavior of vor- 
tices observed in high resolution scanning tunneling mi- 
croscopy images obtained in a slowly decaying magnetic 
field in the range of 0.25-0.75 T^. The images and 
the data derived from them have an unprecedented time 
resolution compared to the average velocity of the vor- 
tices (pm/s). Therefore, the data allows a detailed com- 
parison in terms of track patterns, velocity variations 
and changes in lattice constant. The simulation will 
be mostly used within these field and velocity regimes. 
In this paper we analyze the the general behavior of 
the simulation in terms of collective effects, influences 
of the boundary conditions and defect interactions as 
well as give possible scenarios for the velocity patterns 
observed-. 



The simulation calculates the interaction of an ensem- 
ble of vortices with a static landscape of defects in two 
dimensions. The full vortex- vortex interaction is taken 
into account. The equation of motion for the ith vortex 
is given by: 



Vv - f_^-|)+^FvD (l^i - ^D,j I) + ^global 
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Here Fyy and Fvd are the vortex-vortex and vortex- 
defect interactions, respectively. Fgiobai is an optional 
global driving force. The damping of the vortex motion 
T] is set to 1 (critical damping). 

The main challenge of the calculation lies in the strong, 
long-range vortex- vortex interaction, the large number 
of vortices combined with the desired slow motion of the 
vortex lattice. In a first attempt, the recalculation of the 
vortex-vortex force was controlled for the individual vor- 
tex. The values of Fy.y. were stored in a table. Each 
vortex was allowed to travel a maximum distance dmax or 
evolve for a maximum time tmax — whichever happened 
first — before recalculating Fy.y.. At this time Fy.y. 
is also updated. If the difference to the previous value 
exceeded a threshold, all forces for vortex j were recalcu- 
lated as well. Unfortunately, this approach led to large 
enough accumulative errors to be noticeable in form of 
lattice vibrations (see next section). Thus, the algorithm 
was changed to use a single time control representing the 
minimum value of all vortices. 

Two configurations were commonly used. The first uti- 
lizes periodic boundary conditions and a constant driv- 
ing force. It was mainly used to study the interaction 
of vortices with point defects^'. The size of the sam- 
ple and the number of vortices control the orientation 
of the vortex lattice, especially for the vortex densities 
of interest. There are two basic possibilities for setting 
up a simulation under periodic boundary conditions. Ei- 



FIG. 1: Final states of selected setup scenarios without exter- 
nal driving force, a) periodic boundary conditions, b) periodic 
in X, bound in y, c) repulsive boundaries, d) point defects, e) 
line defects. 



Force 


Function 


/o [nN] 


do [nm] 


Vortex 


Bessel 


30-40 


180 - 200 


Defect 


Gauss 


-40 - 40 


20 - 100 


Boundary 


Exponential 


4000 


50 



TABLE I: Typical parameters used in the simulation. 



ther the number of vortices and the size of the sample 
are matched to form a perfect lattice or dislocations are 
introduced. In general, the number of vortices and the 
sample size were adjusted as to avoid lattice dislocations, 
while maintaining the target magnetic field. This adjust- 
ment leads to a horizontal alignment of the vortex lattice 
(cf Fig. [T^). Furthermore, in the presence of defects, the 
direction of motion tends to align itself with a main axis 
of the lattice. The direction of motion is hence no longer 
identical to the direction of the driving force. 

The second basic setup involved at least two repulsive 
boundaries (cf Fig. [TJd) and c). Here, a disordered layer 
of vortices forms at the edge of the sample making it 
very difficult to avoid lattice dislocations. The lattice 
is driven by removing or inserting vortices at a given 
point which is usually close to the edge of the sample 
and at random time intervals. This procedure generally 
introduces defects that move across the vortex lattice (see 
below). Also, for the field range considered in this paper, 
a significant motion along the sample edge is observed 
owing to the smoothness of the edge potential. 

All forces are modeled as central forces using the form: 
F{d) = fonj^H{d/do). H can be chosen to be a recip- 
rocal, cubic reciprocal, exponential, Gaussian, or Bessel 
function, d is the distance between the two objects, n^* 
its normalized direction. In the case of a line defect or 
a fixed boundary d is the minimum distance, /o and 
d^ are the force constant and decay length, respectively. 
The force and decay constant for the VV interaction are 
calculated from the magnetic field value based on theo- 
retical curves of muon spin resonance measurements on 
NbSes^^ following: 



/o,vv = 



t 



27r/ioA^ 



(2) 



and 



do = A = Ao • ( 1 + 



H 



i^C2(4.2 K)^ 

with <l>o = 2.07 • IQ-^^ Wb, /io = 1-26 • 10"^ Tm/A, 



/ I 



250 pm 



■ 25 pm ■ 



FIG. 2: (Color online) Lattice vibrations caused by numeri- 
cal errors in the simulation. The images show trajectories of 
neighboring vortices. The lateral scale bar refers to the ex- 
cursions. The vortex-vortex distance remains ~ 70 nm. Top 
row: optimized algorithm for maximum evolution times of 20 
s, 2 s, and 0.1 s respectively. Bottom row: common maximum 
evolution time step for all vortices (20 s, 2 s, and 0.1 s). 



Ao = 144 nm, /3 = 1.56 and i^c2(4.2 K) = 2.13 T^. 
t denotes the sample thickness (usually t = 500 /im). 
Typical force function and parameters are summarized 
in table [D 

The setup of the simulation is controlled by an input 
parameter file which allows to easily change the parame- 
ters which include: 

• sample size, number of vortices 

• vortex insertion/extraction rate 

• periodic/fixed boundary condition 

• fixed/time dependent driving force 

• vortex-vortex interaction 

• point/line defects 

• control parameters (<imax, ^maxv) 

Some examples of different simulation setups are dis- 
played in Fig. [H 



III. RESULTS 

A. Collective lattice motion: Thermal noise 

In early versions of the simulation attempts were made 
to optimize the calculation time. The highest computa- 
tion cost lies in the calculation of the vortex-vortex in- 
teraction. By introducing an individual control of the 
evolution time before recalculating the forces on a given 
vortex the performance could be improved by a factor 
of up to ~ 50. However, a problem emerged in form 
of seemingly random lattice vibrations at amplitudes of 
5p ^ 1 nm. In a special set of simulations this motion was 
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systematically studied (cf. Fig. [2]). For large values of 
the maximum time step of tmax = 20 s, amplitudes of up 
to Atyi =0.75 nm and velocities of ^xh < 500 pm/s were 
observed. The values reduced to A^h = 0.25 nm and 
'^Th < 200 pm/s for tmax = 0.1 s. The amplitude and 
velocities were unacceptable, so the code was changed 
to recalculate all forces at a common time interval. For 
large allowed time steps the problem prevailed although 
at smaller amplitudes {A^h =0.75 pm, ^xh ^ 40 pm/s 
and tmax = 20 s) which could be reduced to Axh < 1 am 
and VTh < 40 fm/s at t^ax = 0.1 s. 

The reason for the motion were quasi random kicks to 
the vortex system, which may be due to rounding errors 
accumulated when calculating the sum of forces acting on 
a vortex. This sum has the particular challenge of adding 
numbers of vastly different magnitudes. Under these cir- 
cumstances the order of the summation might lead to 
cases where Fy^y^ ^ Fy^y^ which is likely responsible 
for the observed effect. Fangohr et. al.^^ also showed, 
that under periodic boundary conditions the number of 
repetitions taken into account can have a significant ef- 
fect on the outcome. This simulation originally used only 
the nearest neighbors either in the simulated area or in 
one of the repetitions. The range was later extended to 
the 2nd nearest neighbor which led to a noticeable re- 
duction in the noise level. The later, of cause, depends 
inversely on the size of the simulated area and vortex den- 
sity. It becomes negligible for typical parameters used in 
the simulations described below. 

The motion can be used to estimate an effective tem- 
perature of the system, although the calculations were 
performed at a nominal temperature of K. Since the 
vortices are assumed to be massless the kinetic energy 
cannot be calculated. The variation of the potential en- 
ergy U from its minimum, however, can be expressed as 
an effective temperature: 



1 i:'Li{U{n)-U{r,)) 



N 



(3) 



The potential for a single vortex is usually described by 
a sum over Bessel functions: 



N 



U{r) = -fyryY,Ko 



ry 



(4) 



The potential seen by a single vortex in the direction of 
the nearest and 2"^ nearest neighbor are plotted in Fig. 
[3l The potential around the vortex position is close to 
parabolic. Calculation of fits to a 2"^ order polyno- 
mial (Fig. \^ y = a-\-bx-\-cx'^) confirmed this assumption 
with correlation coefficients of R ^ 0.9997. The coeffi- 
cient c ~ 3.6 • 10^^ dominates the b ^ 3 • 10^ term, and 
since we are interested in the order-of-magnitude of our 
temperature we can rewrite (j3j): 
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FIG. 3: Potential for single vortex in the direction of the 
nearest and 2^^ nearest neighbor at a magnetic field of 0.5 T. 
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FIG. 4: Parabolic fit of the potential near the vortex position. 
The plots are shifted vertically by 0.1 units for clarity. 



Fig. [5] shows the calculated temperatures of the sim- 
ulations shown in Fig. [2] (except for the second one) 
as a function of time. Starting with unrealistic values 
of several thousand K we finally achieved a temperature 
of ~ 100 mK. We used similar control parameters for 
subsequent simulations. Notably, the initial evolution of 
the temperature depends on the initial condition of the 
vortex lattice. This is most easily seen in the bottom 
two curves. For the top curve the initial positions of 
the vortices were chosen with a random offset < 10% of 
the lattice constant while the exact locations were used 
for the bottom curve. In the first case, the lattice cools 
down to the final temperature whereas it initially heats 
up in the second case. In either case it reaches thermal 
equilibrium after about 25 s. 



Parameter 


Nearest neighbor 


2^^ nearest neighbor 


a 


55.7 


55.7 


b 


-1.16 • 10"^° 


6.42 • 10"^ 


c 


3.52 • 10^^ 


3.77 • 10^^ 


R 


0.9998 


0.9996 



^eff 



fyryc 



N 



kB 



N 



TABLE II: Parameters and correlation coefficient of the poly- 
(5) nomial fit to the potential. 
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FIG. 7: (Color online) Local maxima (cf. Fig. [6|) of velocity 
(left) and distortion (right) at different magnetic fields. The 
lines represent fits of an exponentially decaying function to 
the data. 



FIG. 5: Temperature of the vortex lattice as calculated by 
equation \E\ The temperature shows a clear dependence on 
the algorithm and allowed tmax 
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FIG. 8: (Color online) Parameters of fits to an exponentially 



decay (y = yo + Aq - e" 



x/xo\ 



as a function of vortex density 



Distance [um] 



(magnetic field). The lines are a guide to the eye. 



FIG. 6: (Color online) Distance dependence of the velocity 
(top) and lattice distortion (bottom) for a point defect. The 
lines show the envelope and average curves. The dots repre- 
sent local maxima and were used to estimate the decay of the 
variations as a function of distance by fitting an exponential 
decay. 



B. Point defects 

We used single point defects to gain a basic under- 
standing of the vortex-defect interaction. The goals were 
to ascertain the spacial extent of the velocity /distortion 
effects produced and to what degree an arrangement of 
one or more point defects could reproduce the observed 
velocity patterns. For the simulation described here, the 
number of vortices and the size of the area was adjusted 
to achieve a (except for the pinning center) strain and 
defect free vortex lattice, as described above. The den- 
sity was chosen to reflect magnetic fields in the range of 
60-500 mT. The driving force was aligned along one of 
the principal axis of the lattice with a magnitude of 10 
pN, which would lead to a constant velocity of 10 pm/s 
without a defect. 

In the time domain we observed a periodic variation 
of the velocity in the range of 6 - 14 pm/s — largely 



due to the lattice periodicity combined with the periodic 
boundary condition. The distance dependence for a field 
of 125 mT is plotted in figure [H The graph shows the 
average and min/max values of the velocity and lattice 
distortion at a given distance from the defect. To cal- 
culate the width of the affected area we fitted the local 
maxima (black dots) to an exponential decay of the form 
y = yo -\- Aq ' e~^/^° using a least method. The local 
maxima and the fitted curves as a function of the mag- 
netic field are depicted in figure [71 The parameters as 
a function of vortex density are shown in figure [HI The 
lattice shows a clear tendency to soften for lower mag- 
netic fields as the amplitude of the velocity variation and 
the lattice distortion grow. This is not surprising since 
the force constant of the vortex-vortex interaction also 
diminishes with decreasing field according to ([2]). The de- 
cay constant, however, shows only a weak dependence on 
field. The cause for the nonlinearity in xo,v remains un- 
clear at this point, as well as whether the declining trend 
of xo,d continues for lower vortex densities. As noted 
elsewhere^!, although single defects cannot explain the 
time dependent velocity patterns observed in the exper- 
iments, similar lattice distortions were observed which 
indicates the existence of point like pinning centers in 
real samples. 
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FIG. 9: (Color online) Velocity patterns caused by two pe- 
riodic line defects. The image shows the emerging velocity 
patterns. The average velocity of all vortices parallel (top) 
and perpendicular (bottom) to the applied force is plotted 
against time. 



C. Line defects 



Line defects were introduced as a simple representa- 
tion of extended defects in the host material such as step 
edges, dislocations or grain boundaries. The algorithm 
calculates the closest distance of a vortex to the line and 
the resulting force perpendicular to it. Since these types 
of defects are always present in a real material and are 
highly correlated, they present a possible explanation for 
the velocity patterns especially the spikes in the velocity 
we observed^. However, this avenue was not pursued to 
a large extent. The setup using periodic boundary con- 
ditions also requires continuity of the lines at the edge 
and therefore restricts the orientation of the lines. Con- 
sequently, this leads to an overly strong periodic modu- 
lation of the average velocity (Fig. [9)3). 

In the case of repulsive boundaries and moderate force 
values, the line defects have a small influence on the over- 
all pattern (see Fig. [TTl (bottom)). Only strong steps 
show a significant influence on the vortex motion over the 
effect of lattice dislocations (see below). For those cases, 
however, the step increasingly acts as another boundary, 
lining up the vortices and thereby introducing another 
slip plane in addition to the sample edge (cf. Fig. [1^). 
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FIG. 10: (Color online) Enlarged view of a lattice dislocation 
traveling through the field of view from the top right to the 
bottom left corner. The dislocation moves by transposing 
vortices along the direction of travel. 



D. Vortex lattice dislocations 



Dislocations within the vortex lattice proved to be an 
unexpected source of velocity variations. When remov- 
ing vortices to drive the lattice, dislocations are generally 
created near the extraction site or the edge of the sample. 
These dislocation travel through the lattice at moderate 
velocities and lead to speed changes of the underlying 
vortex lattice along the way as shown in Fig. [TOl It 
should be stressed that the dislocations travels much far- 
ther and faster then the individual vortices do. In the 
simulation, the - on a local level - random motion of the 
dislocations lead to seemingly random velocity "spikes" 
of the vortex lattice. For large enough simulated areas 
the lattice should also break up into domains of local or- 
der. This could result in a plate tectonic of the domains 
(getting stuck/unstuck on neighboring domains) which 
in turn could lead to the velocity patterns we observed 
in our measurements. This proposed effect has not yet 
been observed in our simulations. Usually the number 
of vortices is too low to allow '"free floating'" domains 
which are not connected to the sample edge to develop. 

Any defect introduced into one of these or similar large 
scale simulations has only a relatively small influence on 
the overall track patterns that emerge. This is shown 
in Fig. [TTJ after the introduction of two line defects 
the track patterns remain somewhat similar. Due to the 
small area, the overall orientation of the vortex lattice 
is still governed by the boundary. For much larger areas 
(beyond our current computer capacity), this could easily 
change. 
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FIG. 11: (Color online) Lattice dislocation in large vortex 
array (7729 vortices) with repulsive boundaries. Vortices are 
extracted on the right side as indicated leading to a rear- 
rangement of the lattice. The line like structures are vortex 
dislocation traveling through the lattice (see also Fig. [TO]) . 
The lower plot shows the result after adding two line defects. 
Even these correlated defects have little effect on the disloca- 
tion pattern. 



IV. CONCLUSION 



In summary, we developed our own flexible 2D simu- 
lation to calculate the motion of an ensemble of vortices 
under various conditions in a landscape of defects. Based 
upon the average vortex displacement from its potential 
minimum we assert a typical temperature of the simula- 
tion of T ~ 100 mK. We explored the effect of a single 
point defect on the velocity patterns and lattice distor- 
tion and determined the radial extent of the influence. 
The amplitude increases with a decreasing magnetic field 
- in accordance with a softening of the lattice - whereas 
the decay constant only shows a weak dependence. Line 
defects as a model for surface steps or other correlated 
defects were also briefly discussed. 



A new possibility for velocity variations within the vor- 
tex lattice was discovered in form of lattice dislocations 
traveling through the vortex lattice. For larger ensem- 
bles of vortices we could also envision locally ordered do- 
mains moving in unison producing velocity patterns by 
stick-slipping along their neighbors. 
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